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Abstract. We discuss the magnetocaloric properties of gadolinium containing magnetic molecules which 
potentially could be used for sub-Kelvin cooling. We show that a degeneracy of a singlet ground state 
could be advantageous in order to support adiabatic processes to low temperatures and simultaneously 
minimize disturbing dipolar interactions. Since the Hilbert spaces of such spin systems assume very large 
dimensions we evaluate the necessary thermodynamic observables by means of the Finite- Temperature 
Lanczos Method. 

PACS. 75.10.Jm Quantized spin models - 75.40.Mg Numerical simulation studies - 75.50.Xx Molecular 
magnets 



1 Introduction 

The use of magnetic molecules for sub-Kelvin cooling is 
one of the latter ideas in the field of molecular magnetism. 
Originally paramagnetic salts have been employed for this 
purpose which in 1949 led to the Nobel prize for William 
Francis Giauque [J . In recent years several molecular sys- 
tems have been discussed with respect to their magne- 



zos Method (FTLM) [TTIIT^ . In a recent pubhcation we 
demonstrated that this method is indeed capable of evalu- 
ating thermodynamic observables for magnetic molecules 
with an accuracy that is nearly indistinguishable from ex- 
act results [131 . 

The interesting physical question is, which arrange- 
ments of interacting magnetic ions permit a large mag- 
netocaloric effect. For this purpose the low-lying zero- 



tocaloric properties [1[3I1[31[S1[71[H] . Some aspects have turne^gi^i density of states should be high, and an applied field 



out to be of importance if one wants to compete with the 
cooling power of paramagnetic salts: the low-lying energy 
eigenstates should possess large magnetic moments and 
they should be dense, so that one can sweep many of them 
with moderate magnetic fields. The magnetic anisotropy 
should be small, otherwise energy barriers could prevent a 
high density of states. These aspects make gadolinium the 
preferred metal ion: it has a large spin quantum number of 
s = 7/2, a rather small exchange interaction in chemical 
complexes and a negligible single-ion anisotropy, compare 
e.g. Refs. 

From the point of view of theoretical modeling gadolin- 
ium compounds are demanding since the sizes of the re- 
spective Hilbert spaces grow very rapidly with the number 
of involved gadolinium ions due to the large spin quantum 
number. Even when employing all possible symmetries the 
system size is rather restricted since the dimension of the 
largest subspace should not exceed a size of 100,000 in 
order to render a complete numerical diagonalization pos- 
sible [5l[TU]. Fortunately, a very accurate approximation 
has been developed for cases with Hilbert space dimen- 
sions of up to roughly 10^° - the Finite Temperature Lanc- 
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should fan out these levels. Magnetic frustration is known 
to yield dense spectra and to enhance the magnetocaloric 
effect under certain circumstances [T^[T51[T5] . In addition 
frustration, i.e. a situation where not all pairs of spins as- 
sume a coUinear arrangement in the classical ground state 
[17] , may result in singlet ground states which are degener- 
ate. This would have two advantageous consequences: the 
ground state would be connected to an isentrope with S > 
which should allow for very low-temperature cooling. In 
addition a singlet ground state minimizes the dipolar in- 
teractions compared to high-spin molecules where dipolar 
interactions prevent a cooling to very low-temperatures 
|18) . We investigate our hypothesis with the simplest of 
such systems, the tetrahedron. 

It should be mentioned that one can also aim at de- 
signing the intermolecular interactions in order to reduce 
the influence of dipolar interactions [18| . The most simple 
way consists in increasing the mutual distances by means 
of bulky ligands. This route will not be discussed in the 
present article. 

The article is organized as follows. In Section [2] basics 
of the finite-temperature Lanczos method are repeated. 
Section [3] is devoted to the discussion of several recently 
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synthesized gadolinium containing compounds. The arti- 
cle closes with a summary. A technical appendix explains 
a basis coding scheme used in the FTLM. 



2 The finite-temperature Lanczos method 

For the evaluation of thermodynamic properties in the 
canonical ensemble the exact partition function Z depend- 
ing on temperature T and magnetic field B is given by 



z(r,i3) = ^(He-''^ 



(1) 



Here { 1 ) } denotes an orthonormal basis of the respec- 
tive Hilbert space. Following the ideas of Refs. [Tl ] IT2 ] the 
unknown matrix elements are approximated as 



E 

n=l 



V I n(y) )e 



(2) 



For the evaluation of the right hand side of Eq. ^ I ) is 
taken as the initial vector of a Lanczos iteration. This iter- 
ation consists of Nl Lanczos steps, which span a respec- 
tive Krylow space. As common for the Lanczos method 
the Hamiltonian is diagonalized in this Krylow space. This 
yields the Nl Lanczos eigenvectors | n{v) ) as well as the 
associated Lanczos energy eigenvalues Hn'' ■ They are enu- 
merated by n = 1, . . . , Nl- The notation n(y) is chosen to 
remind one that the Lanczos eigenvectors | n{v) ) belong 
to the Krylow space derived from the original state \v). 

The number of Lanczos steps Nl is a parameter of the 
approximation that needs to be large enough to reach the 
extremal energy eigenvalues but should not be too large 
in order not to run into problems of numerical accuracy. 
Nl ~ 100 is a typical and good value. 

In addition, the complete and thus very large sum over 
all states 1 ) is replaced by a summation over a subset of 
R random vectors. These vectors are truly random, they 
do not need to belong to any special basis set. Altogether 
this yields for the partition function 

Z(T,B) « E e-^'S^^'|(nH|.)P . (3) 



R 



u—l n—1 



It will of courseimprove the accuracy if symmetries are 
taken into account as in the following formulation 



jj—l n—1 



(4) 



r labels the irreducible representations of the employed 
symmetry group. The complete Hilbert space is decom- 
posed into mutually orthogonal subspaces 'H{r). 
An observable would then be calculated as 



0{T,B) 



1 



Z{T,B 



E 



dim(H(r)) 
Rr 



Rt 



Nr. 



EE 



x{n{v,r)\o\v,r){v,r\n{v,r)) 



(5) 



This approximation of the observable 0{T,B) may con- 
tain large statistical fluctuations at low temperatures due 
to the randomness of the set of states { | z^, I^)}, but this 
can be cured by assuming a symmetrized version of Eq. ^ 
|19) . For our investigations this is irrelevant. 

In this article the entropy plays a central role, it is 
evaluated as 

S{T, B) = {{H ))/T + kB log(Z(r, B)) , (6) 

with Z{T, B) being calculated according to Q and {{ H )) 

according to ([5|. 

Our very positive experience is that even for large 
problems the number of random starting vectors as well as 
the number of Lanczos steps can be chosen rather small, 
e.g. R « 20, Nl « 100, compare Ref. |T3]. 



3 Magnetocalorics of certain gadolinium 
compounds 

The following spin systems are described by the Heisen- 
berg spin Hamiltonian augmented with a Zeeman term, 



(7) 



Jij is the exchange parameter between spins at sites i and 
j. For the sake of simplicity it is assumed that all spins 



0, this (simple) 



have the same (^-factor. Since 

symmetry is used for the finite-temperature Lanczos cal- 
culations. 

In a process of adiabatic demagnetization the temper- 
ature changes with field according to the following ther- 
modynamic relation: 



dT 

Wb 



T f dS 

c [Wb 



(8) 



Here S denotes entropy, T temperature, B magnetic in- 
duction, and C heat capacity. Besides the heat capacity, 
which in our examples does not vary too much for small 
fields and temperatures T « 1 K, the isothermal entropy 
change (f^)^ has a large impact [8 . At very low tem- 
peratures close to T = the isothermal entropy change 
is of course very large if the ground state of the mag- 
netic molecule possesses a large total spin quantum num- 
ber St, since this corresponds to a theoretical entropy of 
S{T = 0,5 = 0) = kB\og{2St + 1)Q Therefore, naively 

^ Here a short remark concerning the third law of thermo- 
dynamics might be necessary: The third law conjectures that 
at T = the entropy of any system is a universal constant 
that can be taken to be 5 = which in turn means that the 
ground state is non-degenerate. Models, however, can show a 
ground state degeneracy and thus a residual entropy at T = 0. 
In reality tiny interactions might split this degeneracy at their 
energy scale, which means that the residual entropy remains 
at its value down to this scale. 
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thinking, all isentropes with an entropy equal to or smaller 
than this value should run into absolute zero. How close 
they come in reality depends on the very small interac- 
tions that become relevant at very low temperatures. The 
most disturbing interaction for high-spin molecules is the 
dipolar interaction, which also in the case of paramag- 
netic salts limits the achievable temperatures. In the fol- 
lowing we therefore also discuss a possible way out of this 
dilemma: molecules that possess an St = ground state 
with residual entropy due to a ground state degeneracy. 



3.1 Gd4Cu8 & Gd4Ni8 

The M=Cu and M=Ni members of the family of Gd4M8 
molecules were synthesized quite recently [7j. The eigen- 
values of the respective spin Hamiltonians could be de- 
termined numerically exactly for the case of Gd4Cu8, but 
not for Gd4Ni8. In the latter case the Finite- Temperature 
Lanczos Method was employed. 

For Gd4Cu8 the model Hamiltonian ([7| includes the 
following parameters: JcdGd — —0.1 cm"'^, 
>^GdCu = +0.9 cm~^, JcuCu = —8.0 cm~^. The spectro- 
scopic splitting factor was taken as g = 2.0. For Gd4Ni8 
the model parameters were chosen as: JcdGd — —0.1 cm^^, 
JcdNi = -1-0.17 cm~-^, JNiNi = -1-12.0 cm~^. Again we took 
g = 2.0. In the following all other interactions or correc- 
tions such as temperature independent paramagnetism, 
different g factors for different ions or a possible single- 
ion anisotropy in the case of nickel have been neglected. 
Despite these approximations all theoretical curves agree 
nicely with the experimental ones published in Ref. [7]. 




T (K) 



T (K) 



Fig. 1. Theoretical heat capacity per molecule for Gd4Cug 
(l.h.s.) and Gd4Ni8 (r.h.s.) at various magnetic fields. 

Figure [l] displays the theoretical heat capacity per 
molecule for Gd4Cu8 (l.h.s.) and Gd4Ni8 (r.h.s.) at vari- 
ous magnetic fields. The behavior is for both compounds 
qualitatively similar. 

The isothermal magnetic entropy change, compare 
Fig. [2j turns out to be very different; it is much larger 
for Gd4Ni8. The reason is that for Gd4Ni8 the low- lying 
multiplets belong to large total spin quantum numbers 
which leads to larger entropies at low temperatures. This 
is made even clearer in the two following plots display- 
ing the isentropes as function of both temperature and 
magnetic field. 
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Fig. 2. Theoretical isothermal entropy change per molecule 
for Gd4Cu8 (l.h.s.) and Gd4Nig (r.h.s.) for various field 
difrerences:-Z\S'(r, B) = -[^(r, B) - ^(r, 0)]. 
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Fig. 3. Theoretical isentropes for Gd4Cu8. 
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Fig. 4. Theoretical isentropes for Gd4Ni8. 



Gd4Cu8 (Fig. [3]) possesses a non- degenerate St = 
ground state that is separated from a triplet and a quintet, 
whereas Gd4Ni8 (Fig. 4]) has a ground state with St — 22. 
In the latter case all isentropes with S < kg log(45) run 
into absolute zero, which is clearly visible in Fig.[4j On the 
contrary, since Gd4Cu8 possesses a non-degenerate St = 
ground state all isentropes approach temperatures T > 
when B goes to zero. 

Although this behavior suggests that Gd4Ni8 should 
be a very good refrigerant, this does not need to be the 
case. At sub-Kelvin temperatures dipolar interactions be- 
come very important. They prevent a closer approach of 
T = 18\. Dipolar interactions could be tamed by molecules 
that possess an S'f = ground state, but a non-degenerate 
ground state would not be helpful due to its vanishing 
entropy. Therefore, we suggest to investigate molecules 
which have a degenerate - the more the better - ground 
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state with St — 0. A ground state degeneracy can be in- 
duced by frustration, thus a tetrahedron with antiferro- 
magnetic coupling would be a first candidate. 

3.2 A fictitious Gd4 tetrahedron 

About half a dozen Gd4 tetrahedra have been synthesized 
to date, none of them was magnetically characterized [20l 
[ ^ E ^ I ^ IM l l ^ . In the following we therefore discuss the 
magnetic properties of a fictitious Gd4 tetrahedron with 
an exchange interaction of JcdGd = —0.1 cm~^ and g = 
2.0. The magnetic heat capacity (not shown) looks pretty 
similar to those already shown; again C ~ Iks for low 
fields at T w 1 K. 
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Fig. 5. Theoretical isothermal entropy change per molecule 
for a Gd4 tetrahedron for various field differences. 
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Fig. 6. Theoretical isentropes for a Gd4 tetrahedron. The low- 
est isentropes overlap at T = in this linear plot. 
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Fig. 7. Theoretical magnetization contours for a Gd4 tetrahe- 
dron. 



The isothermal magnetic entropy change, compare 
Fig. [5j looks unspectacular. Nevertheless, the isentropes 
shown in Fig. [8] demonstrate an unusual behavior which is 
not obvious at first glance. Since the S* = ground state is 
eightfold degenerate all isentropes with S < /cBlog(8) = 
2.08fcB run into absolute zero. But for a tetrahedron this 
happens in a different way compared to molecules with 
high-spin ground state such as Gd4Ni8. The isentropes of 
high-spin molecules approach zero like a paramagnet with 
a rate of 



dT 



T 
B 



(9) 



i.e. rather steeply. For molecules with a diamagnetic ground 
state zero is approached on a rather "flat" trajectory, e.g. 
from high magnetic fields values like i? = 4TatT = lK, 
compare Fig. [6j 

Figure [7| shows the related magnetization contour plot. 
While decreasing the magnetic field also the magnetiza- 
tion decreases and consequently also the dipolar interac- 
tion. 



so-called three-sublattice antiferromagnet [25]. The non- 
degenerate ground state possesses St = 0, higher-lying 
multiplets are highly degenerate beyond their usual de- 
generacy due to magnetic sublevels. The weak point nev- 
ertheless is the non-degenerate ground state which pre- 
cludes successful cooling, compare isentropes in Fig. [8j 
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Fig. 8. Theoretical isentropes for a Gde octahedron. The low- 
est isentropes overlap at T = in this linear plot. 



3.3 A Gdg octahedron 

As a last example we would like to discuss the octahe- 
dron. This is another interesting structure since it is not 
too complicated to be synthesized and it has an interest- 
ing spectrum. The spectrum of an octahedron is that of a 



In addition to the perfect octahedron we would like to 
discuss a recently synthesized distorted octahedron |27j . 
We assume that this octahedron has an approximate C4 
symmetry. Four spins are arranged on the vertices of a 
square, one at the top, another one at the bottom. A rough 
simulation of the available experimental data yielded the 
following exchange integrals: the interaction between top 
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(or bottom) and every spin of the square Jts — Jhs — 
—0.05 cm~^, between nearest neighbors on the square 
Jss — —0.02 cm""'^, and between top and bottom spins 
Jtb = —0.2 cm^^. The spectroscopic splitting factor was 
taken as g — 2.0. 
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Fig. 9. Experimental and theoretical susceptibility of the dis- 
torted octahedron. Experimental values taken from Ref. [27 . 
The solid curve is the result of a complete diagonalization, the 
dots result from the FTLM. 
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Fig. 11. Theoretical isentropes for the distorted Gdg octahe- 
dron. 

4 Summary and Outlook 

The Finite- Temperature Lanczos Method enabled us to 
evaluate the thermal properties of larger gadolinium con- 
taining magnetic molecules. Due to the large intrinsic spin 
of gadolinium these substances are potentially useful sub- 
Kelvin coolers. A major quality criterion is the achievable 
ground state degeneracy. If such a degeneracy could be re- 
alized for a singlet ground state this would also minimize 
disturbing dipolar interactions. 




B (T) 



Fig. 10. Experimental and theoretical magnetization of 
the distorted octahedron. Experimental values taken from 
Ref. [27]. The solid curve is the result of a complete diago- 
nalization, the dots result from the FTLM. 
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As one can see in Figs. |9] and 10 a coupling scheme 
with the above given exchange interactions yields a very 
good approximation of the experimental data. We can now 
predict how this material would behave as an adiabatic 
cooler. Figure [TT] depicts the isentropes of the distorted 
octahedron. 

Although also this molecule has a non-degenerate dia- 
magnetic ground state the isentropes exhibit a much steeper 
slope, i.e. larger cooling rate, compared to the regular oc- 
tahedron. We conjecture that this results from the fact 
that the non-symmetric, i.e. only C4 symmetric interac- 
tions split the highly degenerate multiplets and thus lead 
to a smeared out density of states. In addition, some of the 
interactions are smaller than in the example or a regular 
octahedron which also reduces the size of low-lying gaps. 



A Basis coding for mixed spin systems 

Since Lanczos iterations consist of matrix vector multipli- 
cations they can be parallelized by openMP directives [55] . 
In our programs this is further accelerated by an analyt- 
ical state coding and an evaluation of matrix elements of 
the Heisenberg Hamiltonian "on the fly" . 

To this end an analytical coding for the product basis 
states 



I mi,. 



, m„, . . . ,TOjv, 



(10) 



is needed. Such a coding was already devised in Ref. [28] 
for equal spins. Here we show that this scheme can be 
easily generalized for spin systems consisting of different 
spins Si, so that — s; < < s^. For encoding purposes, 
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and since m„ can be half-integer, the basis states are usu- 
ally rewritten in terms of quantum numbers — Si — rrii 
instead of m^, where = 0, 1, . . . , 2si. 

The non-trivial technical problem of the coding stems 
from the fact that one wants to use the 5*^ symmetry, 

i.e. work in subspaces ■H(M) of total magnetic quantum 
number M. M assumes values from — MMax up to Mmsx 
with MMax = ^iSi- The basis in the subspace H(M) 
is given by all product states | ai, . . . , a^r ) with M = 
MMax — J2i '^i- ^'-'^ usage in a computer program they need 
to be assigned to integer numbers 1, . . . , dim (H(M)). The 
reason is that one usually does not need the basis only once 
at initialization, but at every Lanczos iteration, since the 
sparse Hamiltonian matrix is not stored, but its non-zero 
matrix elements are evaluated whenever needed using 



{^\H\j] 



(11) 



For a direct coding algorithm of basis states in subspaces 
'H{M) it is advantageous that the sizes of the subspaces 
T-L{M) are known analytically j^^. Thus an array can 
be built at startup that contains for a fixed sequence 
si, S2, . . . , Sat the sizes of these subspaces 7^(M = M^jax — 
A) for given a given number n of spins and A. We will call 
this array D{N,A). It will be used to determine the se- 
quential number of a basis vector in 'H[M). The recursive 
buildup is performed using the following relation between 
the sizes of subspaces 



2s„ 
k=0 



l,A-fc) 



(12) 



If A ^ {0, 1, ... , 2MMax} then D{n, A) = 0. 
For D{n = 1, A = 0, 1, . . . , 2s„) = 1, D{n, A = 0) = 1, and 
D{n,A = l) = n. If A ^ {0, 1, . . . , 2MMax} then 
D{n,A) = 0. 



A.l i 



One coding direction, i ^ | a^^ , . . . , ajy ) , which is the more 
trivial direction, can be realized in several ways. A direct 
algorithm i ^ \a\, . . . ^ a]^ ) using the known dimensions 



of the subspaces 'H{M 
followfl!] 



Ns — A) could be realized as 



m=0 

Ak = A 

do k=N,2,-l 

do n=0,2*s(k) 
if (i.le. (m+D(k-l,Ak-n+l))) then 
BasisVector (k) = n 
Ak = Ak - n 
goto 100 



^ The given code uses FORTRAN notation. Nevertheless, it 
can be easily transformed into C. One should only pay atten- 
tion to the fact that field indices in FORTRAN start at 1 not 
at 0. Therefore, the definition of the second field index of D 
has been modified accordingly. 



100 



else 

m = m + D(k-l,Ak-n+l) 
endif 
enddo 
continue 
enddo 

BasisVector (1) = Ak 



BasisVector contains the N entries a^. This algorithm 
will be made clearer when we explain the inverse algorithm 
below. 
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The inverse direction is actually the nontrivial one, since 
the basis vectors are only a subset of the full basis set 
(10 1. Therefore, for the latter coding direction search al- 
gorithms are often employed, [30] or the two-dimensional 
representation of Lin is used [3T] . 

The position of a basis vector | ai , . . . , a^r ) in the lex- 
icographically ordered list of vectors will be determined 
by evaluating how many vectors lay before this vector. 
For this purpose the known dimensions of the subspaces 
H(M = Mjviax — A) are used again. In a computer pro- 
gram one can evaluate the position i of | ai, . . . , qn ) in 
the list of basis vectors according to 

Ak = A 
i = 1 

do k=N,2,-l 

do n=0,BasisVector(k)-l 
i = i + D(k-l,Ak-n+l) 
enddo 

Ak = Ak - BasisVector (k) 
enddo 

BasisVector contains the N entries a^. If the array of di- 
mension D{N, A) is properly initialized, i.e. the field value 
is zero for non-valid combinations of N and A, then the 
sum can be performed in a computer program without 
paying attention to the restrictions for the indices. 
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